SUBROUTINE DG_ENE
!-----------------------------------------------------------------------
!
! Definition von DG:  Energiegleichung
!

      use physco,   only : z1, z2, z4, z8, z12, pi, sqrt2pi, StBolz, clight, grav
      use primvar,  only : DG, X, XZ, XA, ME, MD, MR, MUr, MUp, MJ, &
                           DG_i, DG_im1, DG_ip1, DG_im2, DG_ip2
      use global,   only : tst, zz, relax
      use config,   only : np, BCflag, rkind
      use advecvar, only : e_adv, de_adv_dim2, de_adv_dim1, de_adv_di, de_adv_dip1
      use geomvar,  only : S_flux, S_vol, S_volA, S_volZ
      use viscvar,  only : divU, divU_flag, alpha_r, q_quad_r, muQ_r, Uterm_r, alpha_p_ene, q_quad_p, muQ_p_ene, Uterm_p
      use matvar,   only : Pgas0Z, dPgas0ZdE, dPgas0ZdD, &
                           TgasZ, dTgasZdE, dTgasZdD, &
                           OPAplkZ, dOPAplkZdE, dOPAplkZdD
      use eddvar,   only : mu_z
      use RBvar,    only : L_star, H_star, f_irr
      use zvar,     only : H_pZ, H_pZ_flag, dH_pZdE, dH_pZdR, H_p, dH_pdE, dH_pdR


      implicit none

      integer :: i


!-----------------------------------------------------------------------
!    Randbedingungen
!-----------------------------------------------------------------------

      ! innere Pseudozellen: i=1 & i=2
!~       if ( X(MUr,3) * tst / X(MR,3) > +1.e-30_rkind ) then  ! Einstroemung ueber Innenrand

         DG(ME,ME,DG_i,1)      =  z1
         DG(ME,ME,DG_i,2)      =  z1

!~       else
!~
!~          DG(ME,ME,DG_i,1)      =  z1
!~          DG(ME,ME,DG_ip1,1)    = -z1
!~          DG(ME,ME,DG_i,2)      =  z1
!~          DG(ME,ME,DG_ip1,2)    = -z1
!~
!~       end if


      ! aeussere Pseudozellen + ueberzaehlige skalare Zelle: i=np, i=np-1, i=np-2
!~       if ( X(MUr,np-2) * tst / X(MR,np-2) < -1.e-30_rkind ) then  ! Einstroemung ueber Aussenrand

         DG(ME,ME,DG_i,np)     =  z1
         DG(ME,ME,DG_i,np-1)   =  z1
         DG(ME,ME,DG_i,np-2)   =  z1

         DG(ME,ME,DG_i,np-3)   =  z1
         DG(ME,ME,DG_im1,np-3) = -z1

!~       else
!~
!~          DG(ME,ME,DG_i,np)     =  z1
!~          DG(ME,ME,DG_im1,np)   = -z1
!~          DG(ME,ME,DG_i,np-1)   =  z1
!~          DG(ME,ME,DG_im1,np-1) = -z1
!~          DG(ME,ME,DG_i,np-2)   =  z1
!~          DG(ME,ME,DG_im1,np-2) = -z1
!~
!~       end if


!-----------------------------------------------------------------------
!    Restlicher Bereich
!-----------------------------------------------------------------------

      do i=3,np-4

DG(ME,MR,DG_im1,i)=-((muQ_p_ene(i)*S_volZ(i)*tst*(z1  - BCflag(i) - BCflag(i-1))**2*(XZ(MUp,i) - &
XZ(MUp,i-1))*z12 *z2 **2*(-(((XZ(MUp,i) + XZ(MUp,i-1))*z12 ) / XZ(MR,i)) + ((XZ(MUp,i) - XZ(MUp,i-1))*z2 ) / (XZ(MR,i+1) - &
XZ(MR,i-1)))*zz ) / (XZ(MR,i+1) - XZ(MR,i-1))**2)
DG(ME,MR,DG_i,i)=-(X(ME,i)*pi*X(MD,i)*(X(MR,i) + X(MR,i+1))) + X(ME,i)*pi*X(MD,i)*(X(MR,i+1) - X(MR,i)) + &
pi*(XZ(MR,i) + XZ(MR,i+1))*tst*z12 *(muQ_r(i)*(z1  - BCflag(i))**2*(-(((XZ(MR,i)*XZ(MUr,i) + XZ(MR,i+1)*XZ(MUr,i+1))) / &
(XZ(MR,i)**2 + XZ(MR,i+1)**2)) + (XZ(MUr,i+1) - XZ(MUr,i)) / (-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 )))**2 + muQ_p_ene(i)*(z1  &
- BCflag(i) - BCflag(i-1))**2*(-(((XZ(MUp,i) + XZ(MUp,i-1))*z12 ) / XZ(MR,i)) + ((XZ(MUp,i) - XZ(MUp,i-1))*z2 ) / (XZ(MR,i+1) - &
XZ(MR,i-1)))**2)*zz  - pi*tst*(XZ(MR,i+1) - XZ(MR,i))*z12 *(muQ_r(i)*(z1  - BCflag(i))**2*(-(((XZ(MR,i)*XZ(MUr,i) + &
XZ(MR,i+1)*XZ(MUr,i+1))) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)) + (XZ(MUr,i+1) - XZ(MUr,i)) / (-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 &
)))**2 + muQ_p_ene(i)*(z1  - BCflag(i) - BCflag(i-1))**2*(-(((XZ(MUp,i) + XZ(MUp,i-1))*z12 ) / XZ(MR,i)) + ((XZ(MUp,i) - &
XZ(MUp,i-1))*z2 ) / (XZ(MR,i+1) - XZ(MR,i-1)))**2)*zz  - XZ(MJ,i)*mu_z*pi**2*(XZ(MR,i) + XZ(MR,i+1))*tst*z8 *zz  + &
XZ(MJ,i)*mu_z*pi**2*tst*(XZ(MR,i+1) - XZ(MR,i))*z8 *zz  - (-(pi*(XA(MR,i) + X(MR,i))) - pi*(X(MR,i) - XA(MR,i)) + &
pi*tst*XZ(MUr,i)*z2 *zz )*e_adv(i) - pi*sqrt2pi*tst*XZ(MUr,i)*z2 *zz *H_pZ(i)*Pgas0Z(i) + pi**2*XZ(MD,i)*(XZ(MR,i) + &
XZ(MR,i+1))*tst*z4 *zz *OPAplkZ(i)*(XZ(MJ,i) - (StBolz*TgasZ(i)**4) / pi) - pi**2*XZ(MD,i)*tst*(XZ(MR,i+1) - XZ(MR,i))*z4 *zz &
*OPAplkZ(i)*(XZ(MJ,i) - (StBolz*TgasZ(i)**4) / pi) + pi*sqrt2pi*tst*(XZ(MR,i+1)*XZ(MUr,i+1) - XZ(MR,i)*XZ(MUr,i))*z2 &
*Pgas0Z(i)*dH_pZdR(i) - f_irr*H_pZ_flag(i)*L_star*tst*((zz *(-(H_star) + H_pZ(i))) / XZ(MR,i)**2 - (dH_pZdR(i)) / XZ(MR,i)) + &
(pi**2*X(MD,i)*XZ(MD,i)*(XZ(MR,i) + XZ(MR,i+1))*tst*(XZ(MR,i+1) - XZ(MR,i))*z4 *(XZ(MJ,i) - (StBolz*TgasZ(i)**4) / &
pi)*dH_pdR(i)*dOPAplkZdD(i)) / (sqrt2pi*H_p(i)**2) - (pi*X(MD,i)*tst*(XZ(MR,i+1)*XZ(MUr,i+1) - XZ(MR,i)*XZ(MUr,i))*z2 &
*H_pZ(i)*dH_pdR(i)*dPgas0ZdD(i)) / H_p(i)**2 - S_volZ(i)*tst*z12 *((muQ_p_ene(i)*(XZ(MUp,i) + XZ(MUp,i-1))*(z1  - BCflag(i) - &
BCflag(i-1))**2*z12 *z2 *(-(((XZ(MUp,i) + XZ(MUp,i-1))*z12 ) / XZ(MR,i)) + ((XZ(MUp,i) - XZ(MUp,i-1))*z2 ) / (XZ(MR,i+1) - &
XZ(MR,i-1)))*zz ) / XZ(MR,i)**2 + muQ_r(i)*(z1  - BCflag(i))**2*(-(((XZ(MR,i)*XZ(MUr,i) + XZ(MR,i+1)*XZ(MUr,i+1))) / (XZ(MR,i)**2 &
+ XZ(MR,i+1)**2)) + (XZ(MUr,i+1) - XZ(MUr,i)) / (-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 )))*z2 *(-((XZ(MUr,i)*zz ) / &
(XZ(MR,i)**2 + XZ(MR,i+1)**2)) + ((XZ(MUr,i+1) - XZ(MUr,i))*zz ) / (-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 ))**2 + &
(XZ(MR,i)*(XZ(MR,i)*XZ(MUr,i) + XZ(MR,i+1)*XZ(MUr,i+1))*z2 *zz ) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)**2) + (z1  - BCflag(i) - &
BCflag(i-1))**2*(-(((XZ(MUp,i) + XZ(MUp,i-1))*z12 ) / XZ(MR,i)) + ((XZ(MUp,i) - XZ(MUp,i-1))*z2 ) / (XZ(MR,i+1) - &
XZ(MR,i-1)))**2*(-(divU_flag(i)*divU(i)*q_quad_p**2*XZ(MD,i)*XZ(MR,i)*z2 *zz ) + &
(divU_flag(i)*pi*q_quad_p**2*XZ(MD,i)*XZ(MR,i)**2*XZ(MUr,i)*z2 *zz ) / (BCflag(i) + S_volZ(i)) + &
(divU_flag(i)*pi*q_quad_p**2*XZ(MD,i)*XZ(MR,i)**2*(XZ(MR,i+1)*XZ(MUr,i+1) - XZ(MR,i)*XZ(MUr,i))*z2 *(-(pi*(XZ(MR,i) + &
XZ(MR,i+1))*zz ) + pi*(XZ(MR,i+1) - XZ(MR,i))*zz )) / (BCflag(i) + S_volZ(i))**2 + &
(alpha_p_ene(i)*XZ(MD,i)*Sqrt(Pgas0Z(i))*dH_pZdR(i)) / Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i))) + (alpha_p_ene(i)*XZ(MD,i)**2*z12 &
*Sqrt(Pgas0Z(i))*dH_pZdR(i)) / (sqrt2pi*(XZ(MD,i) / (sqrt2pi*H_pZ(i)))**1.5*H_pZ(i)) - (alpha_p_ene(i)*X(MD,i)*XZ(MD,i)*z12 &
*H_pZ(i)*dH_pdR(i)*dPgas0ZdD(i)) / (sqrt2pi*H_p(i)**2*Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i)))*Sqrt(Pgas0Z(i)))) + (z1  - &
BCflag(i))**2*(-(((XZ(MR,i)*XZ(MUr,i) + XZ(MR,i+1)*XZ(MUr,i+1))) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)) + (XZ(MUr,i+1) - XZ(MUr,i)) / &
(-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 )))**2*(-(divU_flag(i)*divU(i)*q_quad_r**2*XZ(MD,i)*XZ(MR,i)*z2 *zz ) + &
(divU_flag(i)*pi*q_quad_r**2*XZ(MD,i)*XZ(MR,i)**2*XZ(MUr,i)*z2 *zz ) / (BCflag(i) + S_volZ(i)) + &
(divU_flag(i)*pi*q_quad_r**2*XZ(MD,i)*XZ(MR,i)**2*(XZ(MR,i+1)*XZ(MUr,i+1) - XZ(MR,i)*XZ(MUr,i))*z2 *(-(pi*(XZ(MR,i) + &
XZ(MR,i+1))*zz ) + pi*(XZ(MR,i+1) - XZ(MR,i))*zz )) / (BCflag(i) + S_volZ(i))**2 + &
(alpha_r(i)*XZ(MD,i)*Sqrt(Pgas0Z(i))*dH_pZdR(i)) / Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i))) + (alpha_r(i)*XZ(MD,i)**2*z12 &
*Sqrt(Pgas0Z(i))*dH_pZdR(i)) / (sqrt2pi*(XZ(MD,i) / (sqrt2pi*H_pZ(i)))**1.5*H_pZ(i)) - (alpha_r(i)*X(MD,i)*XZ(MD,i)*z12 &
*H_pZ(i)*dH_pdR(i)*dPgas0ZdD(i)) / (sqrt2pi*H_p(i)**2*Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i)))*Sqrt(Pgas0Z(i))))) - &
(X(MD,i)*XZ(MD,i)*StBolz*S_volZ(i)*tst*z4 **2*OPAplkZ(i)*TgasZ(i)**3*dH_pdR(i)*dTgasZdD(i)) / (sqrt2pi*H_p(i)**2)
DG(ME,MR,DG_ip1,i)=X(ME,i)*pi*X(MD,i)*(X(MR,i) + X(MR,i+1)) + X(ME,i)*pi*X(MD,i)*(X(MR,i+1) - X(MR,i)) - &
pi*(XZ(MR,i) + XZ(MR,i+1))*tst*z12 *(muQ_r(i)*(z1  - BCflag(i))**2*(-(((XZ(MR,i)*XZ(MUr,i) + XZ(MR,i+1)*XZ(MUr,i+1))) / &
(XZ(MR,i)**2 + XZ(MR,i+1)**2)) + (XZ(MUr,i+1) - XZ(MUr,i)) / (-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 )))**2 + muQ_p_ene(i)*(z1  &
- BCflag(i) - BCflag(i-1))**2*(-(((XZ(MUp,i) + XZ(MUp,i-1))*z12 ) / XZ(MR,i)) + ((XZ(MUp,i) - XZ(MUp,i-1))*z2 ) / (XZ(MR,i+1) - &
XZ(MR,i-1)))**2)*zz  - pi*tst*(XZ(MR,i+1) - XZ(MR,i))*z12 *(muQ_r(i)*(z1  - BCflag(i))**2*(-(((XZ(MR,i)*XZ(MUr,i) + &
XZ(MR,i+1)*XZ(MUr,i+1))) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)) + (XZ(MUr,i+1) - XZ(MUr,i)) / (-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 &
)))**2 + muQ_p_ene(i)*(z1  - BCflag(i) - BCflag(i-1))**2*(-(((XZ(MUp,i) + XZ(MUp,i-1))*z12 ) / XZ(MR,i)) + ((XZ(MUp,i) - &
XZ(MUp,i-1))*z2 ) / (XZ(MR,i+1) - XZ(MR,i-1)))**2)*zz  + XZ(MJ,i)*mu_z*pi**2*(XZ(MR,i) + XZ(MR,i+1))*tst*z8 *zz  + &
XZ(MJ,i)*mu_z*pi**2*tst*(XZ(MR,i+1) - XZ(MR,i))*z8 *zz  - S_volZ(i)*tst*z12 *(-((muQ_p_ene(i)*(z1  - BCflag(i) - &
BCflag(i-1))**2*(XZ(MUp,i) - XZ(MUp,i-1))*z2 **2*(-(((XZ(MUp,i) + XZ(MUp,i-1))*z12 ) / XZ(MR,i)) + ((XZ(MUp,i) - XZ(MUp,i-1))*z2 &
) / (XZ(MR,i+1) - XZ(MR,i-1)))*zz ) / (XZ(MR,i+1) - XZ(MR,i-1))**2) + muQ_r(i)*(z1  - BCflag(i))**2*(-(((XZ(MR,i)*XZ(MUr,i) + &
XZ(MR,i+1)*XZ(MUr,i+1))) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)) + (XZ(MUr,i+1) - XZ(MUr,i)) / (-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 &
)))*z2 *(-((XZ(MUr,i+1)*zz ) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)) - ((BCflag(i) + z1 )*(XZ(MUr,i+1) - XZ(MUr,i))*zz ) / (-(XZ(MR,i)) &
+ XZ(MR,i+1)*(BCflag(i) + z1 ))**2 + (XZ(MR,i+1)*(XZ(MR,i)*XZ(MUr,i) + XZ(MR,i+1)*XZ(MUr,i+1))*z2 *zz ) / (XZ(MR,i)**2 + &
XZ(MR,i+1)**2)**2) + (z1  - BCflag(i) - BCflag(i-1))**2*(-(((XZ(MUp,i) + XZ(MUp,i-1))*z12 ) / XZ(MR,i)) + ((XZ(MUp,i) - &
XZ(MUp,i-1))*z2 ) / (XZ(MR,i+1) - XZ(MR,i-1)))**2*(-((divU_flag(i)*pi*q_quad_p**2*XZ(MD,i)*XZ(MR,i)**2*XZ(MUr,i+1)*z2 *zz ) / &
(BCflag(i) + S_volZ(i))) + (divU_flag(i)*pi*q_quad_p**2*XZ(MD,i)*XZ(MR,i)**2*(XZ(MR,i+1)*XZ(MUr,i+1) - XZ(MR,i)*XZ(MUr,i))*z2 &
*(pi*(XZ(MR,i) + XZ(MR,i+1))*zz  + pi*(XZ(MR,i+1) - XZ(MR,i))*zz )) / (BCflag(i) + S_volZ(i))**2) + (z1  - &
BCflag(i))**2*(-(((XZ(MR,i)*XZ(MUr,i) + XZ(MR,i+1)*XZ(MUr,i+1))) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)) + (XZ(MUr,i+1) - XZ(MUr,i)) / &
(-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 )))**2*(-((divU_flag(i)*pi*q_quad_r**2*XZ(MD,i)*XZ(MR,i)**2*XZ(MUr,i+1)*z2 *zz ) / &
(BCflag(i) + S_volZ(i))) + (divU_flag(i)*pi*q_quad_r**2*XZ(MD,i)*XZ(MR,i)**2*(XZ(MR,i+1)*XZ(MUr,i+1) - XZ(MR,i)*XZ(MUr,i))*z2 &
*(pi*(XZ(MR,i) + XZ(MR,i+1))*zz  + pi*(XZ(MR,i+1) - XZ(MR,i))*zz )) / (BCflag(i) + S_volZ(i))**2)) + (-(pi*(XA(MR,i+1) + &
X(MR,i+1))) - pi*(X(MR,i+1) - XA(MR,i+1)) + pi*tst*XZ(MUr,i+1)*z2 *zz )*e_adv(i+1) + pi*sqrt2pi*tst*XZ(MUr,i+1)*z2 *zz &
*H_pZ(i)*Pgas0Z(i) - pi**2*XZ(MD,i)*(XZ(MR,i) + XZ(MR,i+1))*tst*z4 *zz *OPAplkZ(i)*(XZ(MJ,i) - (StBolz*TgasZ(i)**4) / pi) - &
pi**2*XZ(MD,i)*tst*(XZ(MR,i+1) - XZ(MR,i))*z4 *zz *OPAplkZ(i)*(XZ(MJ,i) - (StBolz*TgasZ(i)**4) / pi) - &
f_irr*H_pZ_flag(i)*L_star*tst*(-((zz *(-(H_star) + H_pZ(i+1))) / XZ(MR,i+1)**2) + dH_pZdR(i+1) / XZ(MR,i+1))
DG(ME,MD,DG_im2,i)=-(XZ(ME,i-2)*S_flux(i)*zz *de_adv_dim2(i))
DG(ME,MD,DG_im1,i)=-(XZ(ME,i-1)*S_flux(i)*zz *de_adv_dim1(i)) + XZ(ME,i-1)*S_flux(i+1)*zz *de_adv_dim2(i+1)
DG(ME,MD,DG_i,i)=X(ME,i)*S_vol(i) - pi**2*(XZ(MR,i) + XZ(MR,i+1))*tst*(XZ(MR,i+1) - XZ(MR,i))*z4 *zz &
*OPAplkZ(i)*(XZ(MJ,i) - (StBolz*TgasZ(i)**4) / pi) - (pi**2*XZ(MD,i)*(XZ(MR,i) + XZ(MR,i+1))*tst*(XZ(MR,i+1) - XZ(MR,i))*z4 &
*(XZ(MJ,i) - (StBolz*TgasZ(i)**4) / pi)*dOPAplkZdD(i)) / (sqrt2pi*H_p(i)) + (pi*tst*(XZ(MR,i+1)*XZ(MUr,i+1) - &
XZ(MR,i)*XZ(MUr,i))*z2 *H_pZ(i)*dPgas0ZdD(i)) / H_p(i) - S_volZ(i)*tst*z12 *((z1  - BCflag(i) - BCflag(i-1))**2*(-(((XZ(MUp,i) + &
XZ(MUp,i-1))*z12 ) / XZ(MR,i)) + ((XZ(MUp,i) - XZ(MUp,i-1))*z2 ) / (XZ(MR,i+1) - &
XZ(MR,i-1)))**2*(-(divU_flag(i)*divU(i)*q_quad_p**2*XZ(MR,i)**2*zz ) - (alpha_p_ene(i)*XZ(MD,i)*z12 *zz *Sqrt(Pgas0Z(i))) / &
(sqrt2pi*(XZ(MD,i) / (sqrt2pi*H_pZ(i)))**1.5) + (alpha_p_ene(i)*zz *H_pZ(i)*Sqrt(Pgas0Z(i))) / Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i))) &
+ (alpha_p_ene(i)*XZ(MD,i)*z12 *H_pZ(i)*dPgas0ZdD(i)) / (sqrt2pi*H_p(i)*Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i)))*Sqrt(Pgas0Z(i)))) + &
(z1  - BCflag(i))**2*(-(((XZ(MR,i)*XZ(MUr,i) + XZ(MR,i+1)*XZ(MUr,i+1))) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)) + (XZ(MUr,i+1) - &
XZ(MUr,i)) / (-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 )))**2*(-(divU_flag(i)*divU(i)*q_quad_r**2*XZ(MR,i)**2*zz ) - &
(alpha_r(i)*XZ(MD,i)*z12 *zz *Sqrt(Pgas0Z(i))) / (sqrt2pi*(XZ(MD,i) / (sqrt2pi*H_pZ(i)))**1.5) + (alpha_r(i)*zz &
*H_pZ(i)*Sqrt(Pgas0Z(i))) / Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i))) + (alpha_r(i)*XZ(MD,i)*z12 *H_pZ(i)*dPgas0ZdD(i)) / &
(sqrt2pi*H_p(i)*Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i)))*Sqrt(Pgas0Z(i))))) + (XZ(MD,i)*StBolz*S_volZ(i)*tst*z4 &
**2*OPAplkZ(i)*TgasZ(i)**3*dTgasZdD(i)) / (sqrt2pi*H_p(i)) - XZ(ME,i)*S_flux(i)*zz *de_adv_di(i) + XZ(ME,i)*S_flux(i+1)*zz &
*de_adv_dim1(i+1)
DG(ME,MD,DG_ip1,i)=-(XZ(ME,i+1)*S_flux(i)*zz *de_adv_dip1(i)) + XZ(ME,i+1)*S_flux(i+1)*zz *de_adv_di(i+1)
DG(ME,MD,DG_ip2,i)=XZ(ME,i+2)*S_flux(i+1)*zz *de_adv_dip1(i+1)
DG(ME,ME,DG_im2,i)=-(XZ(MD,i-2)*S_flux(i)*zz *de_adv_dim2(i))
DG(ME,ME,DG_im1,i)=-(XZ(MD,i-1)*S_flux(i)*zz *de_adv_dim1(i)) + XZ(MD,i-1)*S_flux(i+1)*zz *de_adv_dim2(i+1)
DG(ME,ME,DG_i,i)=X(MD,i)*S_vol(i) + (f_irr*H_pZ_flag(i)*L_star*tst*dH_pZdE(i)) / XZ(MR,i) + &
pi*sqrt2pi*tst*(XZ(MR,i+1)*XZ(MUr,i+1) - XZ(MR,i)*XZ(MUr,i))*z2 *Pgas0Z(i)*dH_pZdE(i) - pi**2*XZ(MD,i)*(XZ(MR,i) + &
XZ(MR,i+1))*tst*(XZ(MR,i+1) - XZ(MR,i))*z4 *(XZ(MJ,i) - (StBolz*TgasZ(i)**4) / pi)*(-((X(MD,i)*dOPAplkZdD(i)*dH_pdE(i)) / &
(sqrt2pi*H_p(i)**2)) + dOPAplkZdE(i)) + pi*sqrt2pi*tst*(XZ(MR,i+1)*XZ(MUr,i+1) - XZ(MR,i)*XZ(MUr,i))*z2 &
*H_pZ(i)*(-((X(MD,i)*dPgas0ZdD(i)*dH_pdE(i)) / (sqrt2pi*H_p(i)**2)) + dPgas0ZdE(i)) - S_volZ(i)*tst*z12 *((z1  - BCflag(i) - &
BCflag(i-1))**2*(-(((XZ(MUp,i) + XZ(MUp,i-1))*z12 ) / XZ(MR,i)) + ((XZ(MUp,i) - XZ(MUp,i-1))*z2 ) / (XZ(MR,i+1) - &
XZ(MR,i-1)))**2*((alpha_p_ene(i)*XZ(MD,i)*Sqrt(Pgas0Z(i))*dH_pZdE(i)) / Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i))) + &
(alpha_p_ene(i)*XZ(MD,i)**2*z12 *Sqrt(Pgas0Z(i))*dH_pZdE(i)) / (sqrt2pi*(XZ(MD,i) / (sqrt2pi*H_pZ(i)))**1.5*H_pZ(i)) + &
(alpha_p_ene(i)*XZ(MD,i)*z12 *H_pZ(i)*(-((X(MD,i)*dPgas0ZdD(i)*dH_pdE(i)) / (sqrt2pi*H_p(i)**2)) + dPgas0ZdE(i))) / &
(Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i)))*Sqrt(Pgas0Z(i)))) + (z1  - BCflag(i))**2*(-(((XZ(MR,i)*XZ(MUr,i) + XZ(MR,i+1)*XZ(MUr,i+1))) / &
(XZ(MR,i)**2 + XZ(MR,i+1)**2)) + (XZ(MUr,i+1) - XZ(MUr,i)) / (-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 &
)))**2*((alpha_r(i)*XZ(MD,i)*Sqrt(Pgas0Z(i))*dH_pZdE(i)) / Sqrt(XZ(MD,i) / (sqrt2pi*H_pZ(i))) + (alpha_r(i)*XZ(MD,i)**2*z12 &
*Sqrt(Pgas0Z(i))*dH_pZdE(i)) / (sqrt2pi*(XZ(MD,i) / (sqrt2pi*H_pZ(i)))**1.5*H_pZ(i)) + (alpha_r(i)*XZ(MD,i)*z12 &
*H_pZ(i)*(-((X(MD,i)*dPgas0ZdD(i)*dH_pdE(i)) / (sqrt2pi*H_p(i)**2)) + dPgas0ZdE(i))) / (Sqrt(XZ(MD,i) / &
(sqrt2pi*H_pZ(i)))*Sqrt(Pgas0Z(i))))) + XZ(MD,i)*StBolz*S_volZ(i)*tst*z4 &
**2*OPAplkZ(i)*TgasZ(i)**3*(-((X(MD,i)*dTgasZdD(i)*dH_pdE(i)) / (sqrt2pi*H_p(i)**2)) + dTgasZdE(i)) - XZ(MD,i)*S_flux(i)*zz &
*de_adv_di(i) + XZ(MD,i)*S_flux(i+1)*zz *de_adv_dim1(i+1)
DG(ME,ME,DG_ip1,i)=-((f_irr*H_pZ_flag(i)*L_star*tst*dH_pZdE(i+1)) / XZ(MR,i+1)) - XZ(MD,i+1)*S_flux(i)*zz &
*de_adv_dip1(i) + XZ(MD,i+1)*S_flux(i+1)*zz *de_adv_di(i+1)
DG(ME,ME,DG_ip2,i)=XZ(MD,i+2)*S_flux(i+1)*zz *de_adv_dip1(i+1)
DG(ME,MUr,DG_i,i)=-(S_volZ(i)*tst*z12 *((divU_flag(i)*pi*q_quad_r**2*XZ(MD,i)*XZ(MR,i)**3*(z1  - &
BCflag(i))**2*(-(((XZ(MR,i)*XZ(MUr,i) + XZ(MR,i+1)*XZ(MUr,i+1))) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)) + (XZ(MUr,i+1) - XZ(MUr,i)) / &
(-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 )))**2*z2 *zz ) / (BCflag(i) + S_volZ(i)) + &
(divU_flag(i)*pi*q_quad_p**2*XZ(MD,i)*XZ(MR,i)**3*(z1  - BCflag(i) - BCflag(i-1))**2*z2 *(-(((XZ(MUp,i) + XZ(MUp,i-1))*z12 ) / &
XZ(MR,i)) + ((XZ(MUp,i) - XZ(MUp,i-1))*z2 ) / (XZ(MR,i+1) - XZ(MR,i-1)))**2*zz ) / (BCflag(i) + S_volZ(i)) + muQ_r(i)*(z1  - &
BCflag(i))**2*(-(((XZ(MR,i)*XZ(MUr,i) + XZ(MR,i+1)*XZ(MUr,i+1))) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)) + (XZ(MUr,i+1) - XZ(MUr,i)) / &
(-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 )))*z2 *(-((XZ(MR,i)*zz ) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)) - (zz ) / (-(XZ(MR,i)) + &
XZ(MR,i+1)*(BCflag(i) + z1 ))))) - pi*XZ(MR,i)*tst*z2 *zz *e_adv(i) - pi*XZ(MR,i)*sqrt2pi*tst*z2 *zz *H_pZ(i)*Pgas0Z(i)
DG(ME,MUr,DG_ip1,i)=-(S_volZ(i)*tst*z12 *(-((divU_flag(i)*pi*q_quad_r**2*XZ(MD,i)*XZ(MR,i)**2*XZ(MR,i+1)*(z1  - &
BCflag(i))**2*(-(((XZ(MR,i)*XZ(MUr,i) + XZ(MR,i+1)*XZ(MUr,i+1))) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)) + (XZ(MUr,i+1) - XZ(MUr,i)) / &
(-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 )))**2*z2 *zz ) / (BCflag(i) + S_volZ(i))) - &
(divU_flag(i)*pi*q_quad_p**2*XZ(MD,i)*XZ(MR,i)**2*XZ(MR,i+1)*(z1  - BCflag(i) - BCflag(i-1))**2*z2 *(-(((XZ(MUp,i) + &
XZ(MUp,i-1))*z12 ) / XZ(MR,i)) + ((XZ(MUp,i) - XZ(MUp,i-1))*z2 ) / (XZ(MR,i+1) - XZ(MR,i-1)))**2*zz ) / (BCflag(i) + S_volZ(i)) + &
muQ_r(i)*(z1  - BCflag(i))**2*(-(((XZ(MR,i)*XZ(MUr,i) + XZ(MR,i+1)*XZ(MUr,i+1))) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)) + (XZ(MUr,i+1) &
- XZ(MUr,i)) / (-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 )))*z2 *(-((XZ(MR,i+1)*zz ) / (XZ(MR,i)**2 + XZ(MR,i+1)**2)) + zz  / &
(-(XZ(MR,i)) + XZ(MR,i+1)*(BCflag(i) + z1 ))))) + pi*XZ(MR,i+1)*tst*z2 *zz *e_adv(i+1) + pi*XZ(MR,i+1)*sqrt2pi*tst*z2 *zz &
*H_pZ(i)*Pgas0Z(i)
DG(ME,MUp,DG_im1,i)=-(muQ_p_ene(i)*S_volZ(i)*tst*(z1  - BCflag(i) - BCflag(i-1))**2*z12 *z2 *(-(((XZ(MUp,i) + &
XZ(MUp,i-1))*z12 ) / XZ(MR,i)) + ((XZ(MUp,i) - XZ(MUp,i-1))*z2 ) / (XZ(MR,i+1) - XZ(MR,i-1)))*(-((z12 *zz ) / XZ(MR,i)) - (z2 *zz &
) / (XZ(MR,i+1) - XZ(MR,i-1))))
DG(ME,MUp,DG_i,i)=-(muQ_p_ene(i)*S_volZ(i)*tst*(z1  - BCflag(i) - BCflag(i-1))**2*z12 *z2 *(-(((XZ(MUp,i) + &
XZ(MUp,i-1))*z12 ) / XZ(MR,i)) + ((XZ(MUp,i) - XZ(MUp,i-1))*z2 ) / (XZ(MR,i+1) - XZ(MR,i-1)))*(-((z12 *zz ) / XZ(MR,i)) + (z2 *zz &
) / (XZ(MR,i+1) - XZ(MR,i-1))))
DG(ME,MJ,DG_i,i)=mu_z*pi**2*(XZ(MR,i) + XZ(MR,i+1))*tst*(XZ(MR,i+1) - XZ(MR,i))*z8 *zz  - &
pi**2*XZ(MD,i)*(XZ(MR,i) + XZ(MR,i+1))*tst*(XZ(MR,i+1) - XZ(MR,i))*z4 *zz *OPAplkZ(i)

      end do


END SUBROUTINE DG_ENE
